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I. INTRODUCTION 

During last decade we have witnessed an enormous progress both in low-temperature experimental techniques and in the 
design and better characterization of novel materials and cold atomic systems. These developments allow one to explore the 
quantum world in a more fundamental way. In particular, since interactions between particle constituents can lead to unexpected 
phenomena, one would like to achieve sufficient degree of quantum control to take advantage of it. Theoretical work on strongly 
coupled systems is helping in this regard. For example, research on critical phenomena in quantum phase transitions, where 
the Landau-Ginzburg paradigm of broken-symmetry phase transitions does not apply, shows interesting scenarios. Such is the 
case of the recently proposed deconfined quantum critical points where fractionalized excitations may emerge at criticality with 
observable consequences A crucial theoretical bottleneck, however, is the lack of exactly-solvable interacting many -body 
models, since non-perturbative and non-linear phenomena play a relevant role. 

The main goal of this paper is to introduce a generalization of the Gaudin algebra 01, which we name generalized Gaudin 
algebra (GGA), whose quantum invariants can be exactly diagonalized and may be related to Hamiltonian operators of exactly- 
solvable problems of interacting constituents. By exactly-solvable model we mean a model Hamiltonian whose entire spectral 
problem is reduced to an algebraic one (i.e., it is explicitly diagonalized), a fact that is associated to the existence of a certain 
hidden symmetry in the model under consideration. There are larger classes of Hamiltonians characterized by exact solvability 
of only certain part of their spectra; these are called quasi-exactly solvable yj, and the t-J z chain model is an example |4]. 
Clearly, exactly-solvable models may be used as a starting point to construct many other quasi-exactly solvable models. 

As we will see, we have identified the main operator algebra underlying the integrability and exact solvability of many well- 
known models, thus unifying their description in a single algebraic framework. Simply diagonalizing the quantum invariants of 
the GGA is sufficient to solve all those problems, which include the Bardeen-Cooper-Schrieffer (BCS) |2L Suhl-Matthias-Walker 
(SMW) H, Lipkin-Meshkov-Glick (LMG) 0], generalized Dicke (GD) d, and many others of interest in condensed-matter, 
molecular, atomic and nuclear physics. The basic point is that all these various models, which form the general class of XY Z 
Gaudin models, can be derived using different realizations of the generators of the GGA. For example, the BCS model is 
obtained from the quantum invariants of the GGA after representing their generators in terms of fermionic-pair realizations of 
the generators of su(2). A consequence of this unification is that new exactly-solvable models can be realized after a proper 
representation of the GGA. For instance, one can write down exactly-solvable SU(N) spin and mixed representation models, 
such as spin-fermion, spin-boson or fermion-boson Hamiltonians. 

We start by defining the GGA in section [II] We show how the XYZ Gaudin equation naturally emerges from the Jacobi 
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identity for the generators of the GGA. We also introduce the quantum invariants that will serve as the generating functions for 
all conserved quantities of the generalized (integrable) XYZ Gaudin models. 

In section the XX Z Gaudin equation and the diagonalization of the XX Z Gaudin models are studied. We show a 
family of solutions of the XX Z Gaudin equation, which includes the well-known rational, trigonometric, hyperbolic [Toll 
and the new solution found by Richardson lHUl as especial limits. In particular, we show that the latter can be considered as a 
reparametrization of the other three. The main use of these solutions is to design exactly-solvable Hamiltonians with a large set 
of free parameters, thus providing additional freedom to tune interactions. 

In section HVl we consider two possible realizations of the GGA in terms of the generators of ® j su(2) and su(l, 1) which 
allow us to construct (given the analytic properties of the solutions of the Gaudin equation) general Gaudin model Hamiltonians 
that will be exploited in the rest of the paper. Clearly, from the oscillator realizations of su(2) and su(l, 1) in terms 
of canonical fermions and bosons, one can build several interesting many-body Hamiltonians, including the BCS, SMW, LMG, 
and GD. But one is not limited to these oscillators realizations. Indeed, one can use, for instance, SU (N) or hard-core particles 
realizations to construct new exactly-solvable Hamiltonians fl2il . 

Sections IV1 and [VTI present applications of the algebraic framework to various well-known models. They correspond to 
different realizations of the algebras su(2) or su(l, 1) in terms of canonical fermions or bosons. We start section IVlbv solving 
the BCS pairing models in an arbitrary basis and then focus on the analysis of the BCS Hamiltonian in momentum space. We 
study, particularly, multiband pairing Hamiltonians such as the SMW model which is of relevance for the description of two-gap 
superconductivity in MgB2. In section IVTI we analyze bosonic pairing models of interest in cold atom physics. In particular, 
section |VI CJ concentrates on exactly-solvable two-level boson Hamiltonians, which include the nuclear Interacting Boson Model 
(IBM) 1 13], the two Josephson-coupled Bose-Einstein condensates (BECs) model IU4I1 . and the LMG model. Most importantly, 
we show that the LMG model, used for decades to study phase transitions in finite nuclei, is exactly solvable. 

One would like to understand what are the general differences between mean-field approximations to the XYZ Gaudin 
models, and their exact solution. As we will see, for finite systems the distinction is evident and the character of the solutions 
of the two approaches differs substantially. However, does the difference persist in the thermodynamic limit, i.e. the infinite 
(N — > oo) system-size limit? Expectation values of certain observables (e.g., the BCS gap equation or the occupation numbers) 
will be identical, but other observables may pick up the differences. A question that naturally arises concerns the critical behavior 
of the XYZ Gaudin models. One would like to know, for example, what are their quantum critical exponents. It turns out that 
for certain Gaudin models (e.g., the LMG and SMW model of Eq. i66l ) the critical behavior is mean-field 1151 . This is very 
simple to prove by applying tools from Lie algebras and Catastrophe theory, as developed by R. Gilmore 11611 . The general 
analysis is beyond the scope of the present paper and will be presented in a separate publication. Here, however, we will only 
analyze the quantum phase diagram of the BCS model as a function of the interaction strength and show that the transition 
between the superconducting and Fermi-liquid phases is Kosterlitz-Thouless-like, independently of the space dimensionality of 
the lattice. 

Mixing realizations and representations of the generators of the Gaudin algebra lead to new exactly-solvable models. In 
section lVm we illustrate these ideas by solving three types of many-body models: the GD, an exactly-solvable Kondo-like, and a 
spin-boson models. In this way, one finds the formal algebraic connection between these different physical phenomena and BCS 
superconductivity. Section IVIIII deals with differential operator realizations of the Gaudin generators leading to quasi-exactly 
solvable problems in the continuum. Finally, we show in the Appendix that the weak-coupling limit solutions of the generalized 
XX Z Gaudin models are given by the roots of Laguerre polynomials. 



II. GENERALIZED GAUDIN ALGEBRAS 



A. Commutation relations 



Let us introduce the GGA as the set of operators {S* 

E m 7^ Eg) 



S K (E m )}, with k = x,y, z, satisfying the commutation relations (for 



[S m > S£] — , 

[Sm>SJ = i(y mt S m — X m ( Sf) , 
l s ™. s il = K z mi S m - Y ml Sf ) , 
, Pm; Sf ] = i(X me S m - Z m £ Sf) , 



(1) 



where X ml = X(E m , Eg), Y mt = Y{E m ,E t ), and Z m i = Z{E m , Ei) are antisymmetric (i.e. W(x, y) = -W(y, x)) complex 
functions of two arbitrary complex variables E m , Eg labeled by positive integers m and I, respectively. Equivalently, in terms 
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of the k 



-, z basis (and for E m ^ Eg) 



[S±,Sf]=±2V^(S m 



si), 



(2) 



[SmiS^ ] — ±(V^ £ 



where S± = Sf n ± iS y m , and 



(X„^ ± Y m g)/2. (Notice that Sm^ ' and ^ are non-commuting operators, unless 



X m £ — Y m g.) 

The complex functions X m g, Y m g, and Z m g are taken to have the limiting behavior 

lim eX(x, x + e) — f(x) , lim eY(x, x + e) — g(x) , lim eZ(x, x + e) = h(x) 

£—♦0 e— >0 e^O 



(3) 



where f(x), g(x), and h(x) are non-singular functions. Indeed, X, Y, and Z are complex meromorphic functions having poles 
of order one. In particular, when f (x) = g(x) = h(x) the above commutation relations, Eqs. Q, can be analytically continued 
to the case m = I (i.e., E m — > Ei). For example, 



[S* , S» ] = lim i(Y(E m ,E m + e) S z {E m ) - X(E m> E m + e) S z (E m + e)) = -i f (E m ) 



dS z 

u ~'rr 

dE„ 



Then, 



as! 



f \S K S K 1 = 

tcx cy 1 

Pmi S m ] = — J f(E m ) q£Z 

. [S m i S m ] = — i f(E m ) qjJ^ , 

which together with Eqs. Q form an infinite-dimensional Lie Algebra. 

From the Jacobi identities for the generators of this Lie algebra, for example 

Pn. [S x m , S y e }} + [SI [S«, S« ]] + [S« , [S f y , SfJ] = , 

we obtain, considering the antisymmetry of the functions X, Y, Z, the Gaudin equations fPUl 

ZmlXgn + Z nm Yg n + X nm Y m g = . 

Moreover, the relations (for any pair of indices m, £) 

V" 2 ^2 y2 T 1 

^mf ^ ^mf — 1 1 J A raf ~ 'raf — 1 2 

also result from these identities, where are constants independent of E m and Eg. 



(4) 



(5) 



(6) 



(7) 



B. Quantum Invariants 

Let us introduce the generalized Gaudin field operators 

H(E m ) — H m = S m S m + S/ n S^ + S m S m = — (S m S ro + S m S m + 2 S TO S m ) , (8) 

which act on a carrier space H. These operators are not the Casimir operators of the GGA since they do not commute with its 
generators (E m ^ Eg) 

| [H m , Sf] = ± (V+i {St, S z e } + V me {S£, Sf } - Z ml {Sf,S* m }) , 
[ [H m , S Z {,\ = V m g (Sg S m — S+S^ ) + V mt (S e S m — S m S t ) , 

where {A, B} = AB + BA is the anticommutator. A key property is that these field operators form a commutative family 

[H m ,Hg] = 0, (9) 

therefore, they have a common set of eigenvectors in H and consequently can be considered as a generating function for all 
conserved quantities of quantum integrable systems which will be called generalized XYZ Gaudin models. 
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in. THE XX Z GAUDIN MODELS 

In the following we will concentrate on the diagonalization of the XX Z Gaudin models, i.e. the cases where X m g = Y m i. 
As we will see in the applications, this is the most relevant case from a physics standpoint, and the simpler mathematically since 
the Cartan-Weyl basis is easily defined 12211 . 

A. Solutions of the XX Z Gaudin equation 

The Gaudin equation Eq. reduces to 

ZmiXln + Z nm X(. n + X nm X m £ = 0. (10) 

From this expression, together with the antisymmetry of the functions X and Z, one can derive a parametrization for the 
coefficients Z( n and Xg n : 

v X m iX mn ^ Z mn Z m £ + X mn — Z mn 

~ 7 — ? i "In - — . (11) 

From the latter expression, and Zi n = —Z n i, it follows that 

with r a constant that is independent of any indices. Taking E r as a reference parameter, one can write down 

x r t,x rn r + z rn z r i 

Xln = — _ „ , Zl n = ■ (13) 

Then, the functions X( n and Z( n , which satisfy the Gaudin equations Eq. dl0> . can be written in terms of a limited set of 
parameters s, g, and ti as 



Xi n = g , Z ln =g— — , (14) 

— In t£ ~ tn 

with 

T = sg 2 , U = -g/Z ri , (15) 
where g is a real number and \s\ = or 1. Taking the limit ti — > t„ = x, one finds the limiting behavior defined in Eq. (0 as 

f(x) =g(x) = h(ar) =g(l + sx 2 ). (16) 

In all practical cases one can take the square roots in Eq. (II 4I > to be real and positive (normally any phase can be absorbed in the 
definition of the generators S + and S~). Furthermore, the condition that the resulting exactly-solvable Hamiltonians should be 
Hermitian leads in most cases to the condition that the parameters V and ti be real. In this case the parameter s is either +1,-1 
orO. 

This corresponds to the three cases discussed by Gaudin fl7ll : 



1. Rational: T = 0, s = 0, 



with U = r)i, 
2. Trigonometric: T > 0,s = +1, 



X(rn,r] n ) = Z(m,r] n )=g — - — , (17) 

m - Vn 



X(ve,Vn) = g . i r , Z(r)z,r] n ) = gcot(r]e-r} n ) , (18) 

sm^ - r) n ) 



with ti = tan(^), 
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3. Hyperbolic: T < 0, s = — 1, 

x (vt,Vn) = 9 ■ w 1 r , Z{ve,Vn) =ffcoth(^ - r?„) , (19) 

smh(?7^ - ?7„) 

with = tanh(^). 

Note that for these three parametrizations one finds that the limiting behavior is given by 

f(x)=g(x) = h(x)=g, (20) 
and that the rational model corresponds to the limit rje — ► r) n of both the trigonometric and the hyperbolic model. 



Recently, Richardson has proposed a new family of solutions, given by 



igono 
& 



v , s V 1 + 2aZ t + P z £\/ 1 + 2(XZ n + P z n ry, \ 1 + Ot(Zi + Z n ) + /3ZiZ„ 

X(z e ,z n ) = s-i ,Z(z£,z„) = . (21) 

Zt — Z n Zg — Z n 

Evaluating expression Eq. i ll 2t for this parametrization, one finds that r = j3 — a 2 . Hence depending on the sign of (3 — a 2 , one 
finds that this solution might be expressed as a reparametrization of the rational, trigonometric or hyperbolic models. 
Another useful parametrization is given by 

which can be derived from the hyperbolic parametrization, Eq. dl9> , by taking Cj = e Vi . 

B. Diagonalizing the XX Z Gaudin models 

To define the representation (or carrier) space H of the generalized XX Z GGA we introduce the lowest-weight vector |0), 
such that, 

S-|0) = , S z m |0) = F(E m )\0) \/E m , (23) 

with F(E m ) the lowest- weight function. Thus, the carrier space H is defined as the linear span of the unnormalized vectors 

{|0),S+|0),S+S+|0),-. - ,S+S+---S+|0>,---} . (24) 

We want now to diagonalize the Gaudin field operators. Using Eqs. (0} it turns out that |0) is an eigenstate of H m with 
eigenvalue 

Lo Q (E m ) - F 2 (E m ) - f(E m ) -JL—F(E m ) . (25) 

oE m 

To solve the general eigenvalue problem 

H m \$) =Lj(E m )\<Z>) , (26) 

we propose the Bethe ansatz (M € Z + ) 

M 

i*>=n s ^°)= s ^ s 2 + --- s Mi°)' ( 27 > 



1=1 



and Eq. l l2"6T i is equivalent to 



M 



(H m - w (E m ))\$) = [H m , Yl S+]|0) . (28) 



Thus, the whole problem reduces to compute the commutator 

M M / l-l M 



[H m , II S ^ = E (II S « ) 5 ^ ( ] ' S « ) ( 29 ) 



1=1 l=\ \ n=l 
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whose action upon the state |0) can be written as 



A I 



M I M \ M I M \ 

I ) = E II S r [Hm, St] |0) + - Y. II S r [[ H m, S+] , S+] |0) 



l±n=\ \r(=tl,n) 



which after some algebraic manipulations reduces to 

M I Al 



A I 



E ( r 2Zme F( - E "^ + E z meZ mn |$) + 2J2\ X m i F[E t ) + J2 X ™tZnl \ *tJ°) > 



( = 1 



with *t given by 



*L = (~[[S+(E n ))S+(E m )( [] S+(S n )). 

n=l n=f+l 

Equating to zero all the coefficients in front of ^t m , defines a set of nonlinear coupled equations 



A I 



F(E t ) + Z n e = 

n(^)=l 



£=!,-■■ ,M, 



(30) 



(31) 



(32) 



(33) 



termed Bethe's equations, which determine the set of complex numbers {E m }. Once they are solved, one uses these solutions 
to write down the eigenvalues 



A I 



A I 



w(E m ) = uj (E m ) + E | F ~ 2Zml F ( E ™) + E 

Z m iZ n 



(34) 



e=i 



IV. EXACTLY-SOLVABLE MODELS DERIVED FROM THE GAUDIN ALGEBRA 



Thus far, we have not assumed any special form for the generators of the GGA. Let us consider now a possible realization in 
terms of generators of the su(2) = su(2) © su(2) © • • • © su(2) algebra, which satisfy the relations 



(35) 



with = 5j . The set of indices j will be denoted by the symbol T, whose cardinal is Nq-. Defining the following operators 

in terms of the su{2) generators, 



->m ~~ X m jS^ , S TO — 1 Z m jS? , 



(36) 



one has a possible realization of the generators of the GGA, Eq. Q 1 18]. Notice that the su(2) generators are not constrained to 
be in any particular irreducible representation. Similarly, one can realize the latter in terms of the generators of the Lie algebra 
su(l, 1), homomorphic to su(2), which satisfy 



with (A'j + )^ = K- , obtaining 



s+ = E x ^ K t > s - = - E x ^ K s - s ™ = - V - E 



(37) 



(38) 



jer 



jer 



For the sake of simplicity, we will proceed with the Gaudin operators defined from the (J)i su(2) generators. Extension to 
su(l, 1) is straightforward after application of the non-unitary homomorphic mapping: S 



S+ S~ 



77i m' m 



-S~ S 2 -> S 2 



It is easy to check that the Gaudin field operators, Eq. (|8), are given by 



H m — y ] Z m \Sl + S ' ( Z m iZ m j S?Sj H — -(Sj^Sj + S"; 



ijer 



4 " 



(39) 
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where we assume that X m \ = X(E m , rji) and Z m \ = Z(E m , rji). Bethe's equations, Eq. J33I . are given by 

M 

1 + 2^^+2 ^2z en = 0, i=l,--- t M, (40) 

where dj is the eigenvalue of S?(K?), i.e. S?(K?)\0) = dj|0). In the weakly interacting limit the solutions of these equations 
are given by the roots of associated Laguerre polynomials (see Appendix [A}. This establishes a one-to-one correspondence 
between the eigenstates of the non-interacting and the weakly-interacting models, which proves that the Bethe ansatz covers all 
eigenstates and does not contain any spurious solutions for finite values of g. 

In the strong interaction limit, g — > ±00, some of the variables diverge to infinity, where again they can be related to roots of 
the associated Laguerre polynomials, scaled with a factor q. However, some of the roots can remain finite. These finite roots are 
equivalent to the solutions for the Gaudin spin magnets 11711 and can be related to the elementary excitations of the BCS-model 
in the canonical ensemble 12611 . 

From the analytic properties of the X and Z function matrices (771 ^ r/j) 

I ^X mi X mi = ffaJJTy , / ^Z mi Z mj = ffa,);j u , / ^Z mi = ffa) , (41) 



where T{ is a contour in the complex-i? m plane encircling 771 . In this way, one can write down the constants of motion R\ = 




for t sit (2) and ® t su(l, 1), respectively, with eigenvalues r\ = jp_ ^foj{E rn ) 

ri = di I 1 + Zu + 2 d > Zi > ■ (44) 

A class of Gaudin model Hamiltonians can be written as Hq = ^\ E{Ri, i.e. 

H G = £ £i <?f + JL £ (x lJ (8+8f+Sr8+)+2Z^Sf) , (45) 

Hq = JL J2 (x ii (K+Kr+K.-K+)-2Z ii K?K?) , (46) 

for ®j su(2) and ®j sit(l, 1), respectively. In the equations above e\ is an arbitrary real number, X-^ = (gj — Ej)Xij/g, 
Z;j = (si — Ej)Zij/g are real symmetric matrix functions (ATy = X(r]i, rjj), Z\- } = Z(r/i,r]j)), and g = gN is a c-number of 
order of magnitude unity because of thermodynamic stability reasons. Notice that, since Ei and rji are, in principle, independent 
parameters, one may take advantage of this freedom to write down different kinds of mode-dependent interactions (see section 
IVIB> . Moreover, e\ and rji must be chosen real for Hq to be Hermitian. Clearly, there are other classes of Gaudin models that 
involve higher-order combinations of the integrals of motion Ri. 

Linear combinations of the i?i's will be used in the next sections to derive different exactly-solvable model Hamiltonians. 
Different models result from using different realizations of su(2) (or su(l, 1)). In the following we will use (canonical) fermion 
and boson realizations, though, one could have used many others Jl2il . such as SU(N) spins or hard-core particles, leading to 
new exactly-solvable problems all of them having the same algebraic root. For example, for SU (2) in the spin S = 1 irreducible 
representation one can write down Eq. ( 145 \ (for the rational case) in terms of sit(3) generators {[i^v — 0,1, 2) in the 
fundamental representation as 11211 

h g = aisl 1 - + J y isr s r - srs^) , (47) 

i ij(i^j) 

where Jy = (ej — £j)Xy/2, and S uft are the generators of su(3) in the conjugate representation. 
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V. 0; su(2) FERMIONIC REPRESENTATION MODELS 
A. BCS-like Models 

Not many models in condensed matter physics have attracted that much attention as the Bardeen-Cooper-Schrieffer (BCS) 
model of superconductivity @|, a remarkable phenomenon discovered in 1911 by Gilles Hoist and Kamerlingh Onnes fiUl . and 
which is characterized by vanishing electrical resistance and perfect diamagnetism ll20ll . Soon after the introduction of the BCS 
model in condensed matter, Bohr, Mottelson and Pines 1I2H1 applied the BCS theory to the description of pairing correlations in 
finite nuclei. The BCS or pairing Hamiltonian is given by 

Hbcs = £ i + E 9iv' A^^vo- ■ (48) 

1 krl'tr' 

The operator c\ a (c 1(T ) creates (destroys) a fermion in the state ler, where a is the third projection of the internal spin degree 
of freedom S, 1 refers to all other quantum numbers needed to specify completely the state, and n\ = Cj & is a number 
operator. Though, in principle, the state ler could be an arbitrary conjugate state to ler (only a bijective relation between conjugate 
pairs is required), we will restrict here to time-reversal conjugate pairs. Under time-reversal (effected by an antiunitary operator) 
the position operator stays unchanged while the (linear or angular) momentum or spin operators change sign. Thus, the time- 
reversal transformation of single-particle states is specific to the choice of basis. For example, the eigenstates of a generic angular 
momentum operator J, labeled as \jm), transform as 

\J^) = (-l) j - m \j-m). (49) 

Similarly, the time -reversal transformation of an annihilation operator in a basis of spin and linear momentum is cj— = 
(— l) s ~ CT c_k-<r, while for a basis of spin and position = Cro = (— l) s_(7 c r _ <7 -. For the sake of clarity, we will assume 
a position basis in the following such that the time-reversal operation will be referred exclusively to the internal spin part of the 
states, i.e. = c^. 

The pairing Hamiltonian ( |48> with uniform couplings gfj? = g/A has been solved exactly in full generality by Richardson in a 
series of papers in the sixties 12311 . This important development escaped the attention of the condensed matter and nuclear physics 
communities until very recently, when the Richardson's works were rediscovered in the study of ultrasmall superconducting 
grains. In order to regain the exact solution we will now present a specific representation of the su(2) generators in terms of 
fermions 

a a 

where the operator Tj + creates a pair of fermions in time-reversal states and £!i = 2t\ + 1 is the degeneracy of the state 1 related 
to the pseudo-spin of the state 1. It can be readily verified that the three operators {rj^, rf } satisfy the su(2) algebra i35\ . 

The integrability of the BCS Hamiltonian d48> was recently demonstrated ll24ll. It was shown that -Hbcs can be written as a 
linear combination of the integrals of motion of the rational family with = = g/(£i — £j) < l 1 71 



# bcs = ei Ri + c = E ei ( 2r i z + 5 n o + 9 E t i +t i 7 • 



(51) 



1 1 if 

The complete set of eigenstates of the pairing Hamiltonian are given by the product wavefunction 

M 

S+ \,A . S+ = V X™, r+ = V 

Em, — 2ei 



i*> = n s ™ i"> > s ™ = e x ^ ^ - e r . < 52 > 

m=l 1 1 

where \v) = \v\,v 2 • ■ • , vl) , with L the total number of single particle states, is a state of v unpaired fermions (v = J2i ^1) 
defined by 

rf \v) = , 711 1^) = v\ \v) ■ 

The quantum numbers v\ are often referred to as Seniority quantum numbers in the nuclear physics literature. 

The total number of particles is N = 2M + v, with M the number of Cooper pairs. Each eigenstate d52t is completely defined 
by a set of M spectral parameters (pair energies) E m which are a particular solution of the Richardson's equations 

l + iy ni - 2 " 1 +2g Y 1 =0. (53) 
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The eigenvalues of the BCS Hamiltonian are 

L M 




E=^2siPi+J2 E m- (54) 

1=1 m=l 

One can easily relate the spectra of the repulsive (g > 0) and attractive (g < 0) cases: If one performs the following canonical 
particle-hole transformation 

I — ^ 

(55) 

which is not a symmetry (although the interaction term is invariant), -Hbcs (.<?) transforms as 

ffBCs(s) E " l£l - H bcs(-9) , (56) 
l 

indicating the relation between the two spectra. 

In recent years, the exact solution of the BCS Hamiltonian has been recovered in the study of ultrasmall superconducting (for 
a review see |3). The specific Hamiltonian for grains assumes a set of L equally spaced doubly-degenerate single particle 
states. Implying that fl\ = 2 and E\ = 1, with 1 = 1, 2, • • • , L. The Richardson's equations d53l > reduce to 



L M 

•^E^— -t+%> E tf—w.=°> (57) 



1=1 t&m)=l 

with i/\ = 0, 1. The ground state for an even number of particles N is in the sector of no broken pairs, v\ = for all 1, while for 
odd JV, V\ = 1 for 1 = (N + 1) /2 and zero otherwise. In other words, the Fermi level is blocked by a single particle, excluding 
it from the active space as can be seen from the second term in equation i57\ . The additional gap at the Fermi energy due to the 
blocking of this level is at the origin of the odd-even difference observed in the tunneling spectra of sm all g rains. The excited 
states of the model are either collective states (pairing vibrations) within the same Seniority subspace II26I1 . or non-collective 
broken pairs l27ll . 

While this pairing model has found great success describing the physics of ultrasmall grains, it can be likewise applied to 
axially deformed nuclei with non equally spaced single particle levels. The reason that prevented its use in standard nuclear 
structure calculations for so many years, was the lack of an efficient numerical procedure to solve the equations d57i for a large 
number of non-equally spaced levels. While in the equally-spaced case the method proposed by Richardson 12811 allowed the 
treatment of systems with ~ 10 3 particles ll29ll . the singularities arising in the numerical solutions of the equations with non- 
equally space levels are difficult to treat. Recently, it has been proposed a new numerical procedure to avoid the singularities 
which seems to be very promising ll3(ill . and it might open the scope for applications to several quantum systems. 

The original BCS model for superconductivity ja was introduced in the context of bulk metallic superconductors. In this case, 
electrons (spin- 1/2 fermions) are confined in an arbitrary dimensional box with periodic boundary conditions with single-particle 
states of the Bloch type. Pairing occurs in momentum k-space. For pairing in the singlet s -wave channel and Cooper pairs with 
zero momentum (k f , k J.) the BCS Hamiltonian can be written as (n^ a = cj^c^ with a =|, J.) 

■Hbcs = E £kCT nkCT + -9 E c kT c -ki c -k'j. c k'T 

k<r kk' 

= E[ £k (" k T + n -kl) -3"kT"-kJ +fE c kT c -ki c -k'i c k'T > (58) 

k kk' 

where the prime in the first double sum means that the terms k = k' are omitted, and cS creates an electron with momentum k 
and spin a. It has been assumed time-reversal invariance, i.e. £kt = £-kj. = £k- The relevant su(2) algebra in this case is 



1 



kr _ ki - (T k ) j Tk = ^("kT + »-ki - 1) . (59) 

where ilk = 2, i.e. the single particle states k | and — k j are degenerate. However, e k may, in principle, differ from e_ k - The 
Hamiltonians of Eqs. (1581 and ( 15 1 1 are dynamically equivalent. To see this let us rewrite Eq. ( 1581 in terms of the pseudospin 
operators r 

^BCS = Et £ k(2^ + 1) - g{ri + 2K) 2 )] + g^ r+r^, . (60) 

k kk' 
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FIG. 1: Quantum phase diagram of the BCS model as a function of the interaction strength g in the thermodynamic, N — > oo, limit. The 
insets represent the single-particle occupation number rik in each quantum phase. Notice that, for positive g, the Fermi liquid quasiparticle 
renormalization factor Z* is unity regardless of the magnitude of the interaction. 

It can be easily shown that the operators ^ k Tjjf, and X)k( T k) 2 WQ conserved quantities, i.e., commute with T/bcs- Thus, up to 
an irrelevant global constant, 

H BCS = £k ( 2T k + !) + 9 T k T k" . ( 61 ) 

k kk' 

which is clearly equivalent to Eq. ( 15 U . The eigenvalues of the BCS Hamiltonian, Eq. ( I58> , are given by Eq. i54i where the 
parameters E m are the solutions of the Richardson's equations, Eq. d53l >. with ilk = 2. One needs to take into account the fact 
that for each k there is a k in those sums. Moreover, if the crystal has space-inversion symmetry £k = £-k- This additional 
symmetry, which converts each single-particle level into a four-fold degenerate one, may have dramatic consequences. For 
example, the numerical solution for the ground state in the BCS case is free of singularities due two the fact that the pair energies 
E m come in complex conjugate pairs for any value of the coupling strength g. 

In previous work 1121 . a gauge ST/ (2) symmetry was identified. The SU(2) symmetry generators are the local operators 

S k = c icT c -ki = ( S k ) f . = \( n ^1 - ™-kl) , (62) 

which commute with the pseudo-spins Tk, i.e [S£, t£,] = 0, for /Lt, v = ±, z. This symmetry amounts to the conservation of 
the charge parity per mode pair (k f, — k J.). Indeed, this local symmetry is responsible for the Pauli blocking of the (unpair) 
singly-occupied states. We would like to emphasize that all the symmetry analysis applied to Eq. fl58l > is also applicable, after 
proper rewritting of the symmetry operators, to Eq. ( I51l l. 

It is interesting to analyze the quantum phase diagram of the BCS Hamiltonian TZbcs as a function of the coupling strength 
g = gN. To this end one needs to study the behavior of the quantum correlations of the ground state in the thermodynamic 
limit. It has been shown, under quite general assumptions, that the Bethe equations of the integrable BCS Hamiltonian in the 
thermodynamic limit are the BCS equations l3lll32l . The condensation energy for attractive pairing in this limit is E con ^ = 
— =^P-e ' 9 where d is the mean level spacing (~ Vol" 1 ) and u; d is the Debye frequency cutoff. For repulsive pairing E con ^ = 0. 
It turns out that there is a quantum phase transition between a BCS superconductor (broken U(l) symmetry) and a Fermi liquid 
of a peculiar type at g = (see Fig. It is important to emphasize that the ground state energy has an essential singularity 
at g = CP, implying that it is an infinite-order (Kosterlitz-Thouless-like) quantum phase transition but with a broken U(l) 
symmetry. (Notice that this result is independent of the space-dimensionality of the problem.) We have numerically solved the 
Bethe equations and found that the Fermi liquid has quasiparticle renormalization factor Z* = 1 independently of the magnitude 
of g; moreover, it displays enhanced superconducting fluctuations but it is not a superconductor. The fact that Z* = 1 has been 
previously remarked in 13311 using functional integrals. 
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The exactly solvable Hamiltonian Eq. J58I can be generalized to the case of multibands in the following way 

#BCSn = XI £ k KkT + ™n-kl) + ^ C nkT C Lkl C n'-k'l C n'k'T » ( 63 > 



nkn'k' 



where n represents the band index, n n k<r = c nkcr c nko-' an( ^ Snn' = ffn'n- The prime in the sum means that the terms (n, k) = 
(n', k') are excluded. Global symmetries of the model include r z = J2 n k r nk anc ^ Xm k( r nk) 2 w ^ tn r nk = 5( n nkT +^n-kj. — 1)- 
The local SU(2) symmetry has as generators 

S nk = 4kT C n-k| = C^nk) 1 > ^nk = T^nkT - ™n-kj.) • (64) 

Clearly, the BCS Hamiltonian of Eq. J58fc is a particular case of Eq. (|63} for a single band (n = 1). It is straightforward to see 
that i?BcSn is exactly solvable for interactions 
SW = <5nn'5n (decoupled BCS bands) 
2.- .gnn' = 5 (effective one-band BCS model) 

For the particular case of two bands, the Hamiltonian Eq. (1631 might be of interest to describe the phenomenon of two-gap 
superconductivity recently observed in materials like MgB2. We may consider here an SU (2)i <g> SU (2)2 structure, with each 
SU(2) n generated by the elements 



T " - 4kT c n-k4 - ( T n ) ' T n - ''r 



The case of two fiat bands e = — = e\ with equal diagonal interaction terms, i.e. g n = g 2 2, can be easily shown to be 
exactly solvable: By using the quantum invariants of the XX Z Richardson-Gaudin (RG) models 

R n =t* + 2J2 [^f-{r+T- f + T-T+) + Z nn ,T*T* n ,\ , (65) 

it can be shown that the two-band pairing Hamiltonian is equivalent (up to an overall constant) to 

#BCS2 = 2e (r| - rf) + g n (r+rf + t+t^) + 912 (r 2 + rf + t+t^) , (66) 

where <7n = 4eZ2i and 512 = ^eX^i are two arbitrary real numbers (with the parametrization of Eq. dl4l> . g\2~9ix = (4£ff) 2 s)- 
To arrive to expression J66l > we have used the Casimir invariants 

\{<r- + r-r+) + (r 2 ) 2 = S(S + 1) (67) 

together with the conservation of r 2 = rf + rf and (r* ) 2 + (r| ) 2 . This is the Hamiltonian originally proposed by Suhl, Matthias 
and Walker |(J as an extension of the BCS model, to include situations where the scattering between electrons from different 
bands contributes substantially to the resistivity in the normal state. 

The (unnormalized) eigenstates of Hbcsi are given by (with t\ = —r\ and t<i = rf) 

M ( 1 1 \ 

where the spectral parameters Ee satisfy Bethe's equations (d± = di ± cfe, and 2d 1 ( 2 ) = ^1(2) — ^1(2) /2) 

£ + v9l2d+ E t -2e 9 d (l + sEt) +2g£ ^ l+^ =0 , m 
A? - £ i E n 

The corresponding eigenvalues can be easily obtained from those of the integrals of motion. 

More complex situations arise in the application of the BCS model to the spherical nuclear shell model or to finite lattices. 
As an example in nuclear physics we will consider the semi-magic Sn isotopes . These series of nuclei can be modelled by a set 
of valence neutrons occupying the single-particle orbits in the N = 50 — 82 shell interacting with a residual BCS Hamiltonian. 
In Table|l]we show the experimental single particle energies and the corresponding degeneracies in the spherical single particle 
basis. 
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TABLE I: Single particle energies and degeneracies for the Tin isotopes in the N = 50 — 82 shell. 



s. p. level 


d 5 /2 


g7/2 


Sl/2 


d 3 /2 


hn/2 


s. p. energy (Mev) 


0.0 


0.22 


1.90 


2.20 


2.80 


s. p. degeneracy 


6 


8 


2 


4 


12 



Richardson's equations for this case have non-equally spaced levels and variable degeneracies. While the solution to this 
problem has been found using standard techniques l34ll . larger systems would require more sophisticated methods l30ll to avoid 
the singularities. 

A quite similar situation arises in solving BCS Hamiltonian in finite two dimensional lattices of size L x L |'9|]. The single 
particle energies in units of the hopping matrix element are Sk = — 2 (cos k x + cos k y ), with k p = 2-Kn p / L and —L/2 < n p < 
L/2. Table im shows the single particle energies and degeneracies for a 6 x 6 lattice. 

Numerical applications of the RG models to fermionic problems in nuclear physics and condensed matter have been concen- 
trated on the BCS Hamiltonian with uniform couplings. The use of the hyperbolic model, a particular solution of the XX Z 
generalized Gaudin models, with non-uniform coupling strength has been suggested in Ref. l35ll to describe the physics of 
multigrain systems, though no practical applications has been carried out so far. 



TABLE II: Single particle energies and degeneracies for electrons in a 6 x 6 square lattice. 



s. p. energy 


-4 


-3 


-2 


-1 





1 


2 


3 


4 


s. p. degeneracy 


2 


8 


8 


8 


20 


8 


8 


8 


2 



Another exactly-solvable model with a separable pairing interaction (SPI) was proposed by Pan, Draayer and Ormand I3f3l . 
The Hamiltonian has degenerate single-particle energies but some structure in the pairing interaction 

H S Pi=eJ2ni+ J2 Sjj'^m^^^jw^'™' ' (70) 

j jj',m,m' 

with <7jj/ = gcjc'- y This model can be derived from the model of Eq. ( 145 \ using the parameterization of Eq. Mil and taking 
= c? . Inserting the number operator divided by the number of particles and adjusting the interaction strength, one can cancel 

out the one-body and two-body diagonal parts in the Hamiltonian. 

Apart from its relevance in the nuclear shell model, the SPI model has also been used in connection with atomic BECs l37ll 

and for establishing variational lower bounds on the energy of general two-body Hamiltonians l38ll . 

B. Particle-hole-like Models 

It is clear that what is behind the exact solvability of these different models is a GGA and the existence of certain quantum 
invariants. Different representations of the Gaudin operators lead to different models but all of them with the same dynamics. In 
this section, we continue with su(2) fermionic representation models. 

In previous section we have written down BCS-like models using an su(2) representation in terms of pseudo-spins r. We have 
also seen that there is another fermionic representation for su(2) in terms of the generators of Eq. \62\ . The natural question 
that arises is: Can we write down sensible exactly solvable models of interacting fermions in terms of this su(2) representation? 
and the simple answer is yes. 

The Lipkin-Meshkov-Glick (LMG) model Q was originally introduced to study phase transitions in finite nuclei. The model 
considers N fermions distributed in two Af-fold degenerate levels (termed upper and lower shells). The latter are separated by 
an energy gap e 

k<r kk'o" kk'o" 

with the quantum number a = ± labelling the level. In Eq. ( 17 H the interaction term V scatters a pair of particles across the 
Fermi level, i.e. it is a two particle-hole interaction, while the term W exchange particles in the two levels. Upon introducing 
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the collective particle-hole operators 

S+ = E 4+Ck- = (S~) f = J E CTC k^ kff , (72) 



2 

k kff 

which satisfy the su(2) commutation relations, Eq. (17 1 i may be rewritten as 

H LMG = sS z + ^(S + S+ + S-S~) + ^(S + S- + SS+) . (73) 

As defined by Eq. I73i . H\_mg is invariant under the inversion symmetry operation / that transforms (S x ,S y ,S z ) >— > 
(S X ,-S V ,S Z ), and it also commutes with the (Casimir) operator S 2 = (S + S~ + S-S+)/2 + (S z ) 2 . Thus, Eq. {73) is 
equivalent to 

H lMG = eS z + ^(S + S + + SS~) - ^S Z S Z + ^S 2 . (74) 

The Hamiltonian -Hlmg nas a band matrix representation in an su(2) basis, and it can be easily diagonalized for large values of 
N. As such, the model has been used as a testing ground for many-body approximations in nuclear physics. More recently, the 
simplicity of the model and the fact that it can be interpreted as a Heisenberg chain with long range exchange interactions, made 
it fashionable to study relations between entanglement and quantum phase transitions. 

We will show in section fVI C 2l that the LMG model is exactly solvable. But before consider the modified problem 



H 9-b = E Y <Tt hr°k* + 7^7 E C L C l-o C k'<r C k-a , (75) 



£k t W_ 
k(T kk'cr 

where a = ± may be now interpreted as a band index. Let us introduce the following commuting su(2) algebras 



= 4 + c k _ = (5 k )t , = l(n k+ - n k _) , (76) 
t£ = 4+4- = far ) f , < = i(n k+ + n k _ - 1) , (77) 



2 

in terms of which i? p -h can be written (up to an irrelevant constant) as 



ffp-h = E £ ^k + f E ■ (78) 



AT 

k kk' 



lf p _h is dynamically equivalent to -Hbcs an d, thus, it is also exactly solvable. 



VI. 0, su(l, 1) BOSONIC REPRESENTATION MODELS 
A. Bosonic BCS-like Models 

The boson BCS or pairing Hamiltonian can be written in complete analogy as in the fermion case, Eq. ( I48> , as 

#BBCS = E £l ni + f E b H V & 1' ' ( ?9 > 
1 11' 

where b\ (by) creates (destroys) a boson in the state 1, and n\ = b\b 1 is the number operator. For simplicity, we will consider 
here scalar bosons, but an arbitrary internal spin and be easily taken into account as in the case of fermions. The label 1 is 
a short-hand notation for a set of quantum numbers; for example, the states of a 3D isotropic harmonic oscillator potential 
are labeled by 1 = (nlm), where n is the oscillator quantum number, / is the orbital angular momentum and m its third 
projection. In Eq. ( I79> T refers to the time-reversed state of 1. Following j49t the time-reversed annihilation boson operator is 

= ^rdrn = (~ m ^nl-m- 

In the following we will be concerned with spin scalar bosons, a possibility that cannot be realized in fermionic models. It is 
worth emphasizing, however, that the exact solution for boson systems can easily incorporate the spin degree of freedom (integer 
spin), and there might be important applications for spinor BECs l39ll not explored so far. 
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Once again, Richardson 14011 determined the complete spectrum of the boson BCS Hamiltonian of Eq. &791 . This work also 
escaped the attention of the physics community until very recently, when the model was shown to be quantum integrable |9]. 
Exactly-solva ble gene ralizations of the uniform pairing Hamiltonian were proposed and subsequently applied to various finite 
Bose systems |4ll|4j. In analogy with the fermionic systems presentation of previous sections, we will first introduce a specific 
representation of the su(l, 1) generators 

K = \^ = {K^ , K{ = 1^ + ^01, (80) 

where the operator Kf creates a pair of fermions in time-reversal states and Q\ = 2K\ + 1 is the degeneracy of the state 1 related 
to the pseudo-spin K\. The operators of Eq. (I80i satisfy the su(l, 1) commutation relations 

The BCS Hamiltonian -Hbbcs can be derived from the rational family dl7> as a linear combination of the integrals of motion 

#bbcs = E ei Ri ( £ ) + c = E £ ^ 2K t - + a E K ? K v ■ (81) 

1 1 11' 

The complete set of eigenstates of this model are given by the product wavefunction 

M 

I*) = J! S » I") . S ™ = E X «* K \ = E O g _ F K l , (82) 

m=l 1 1 

where \v) = \v\, v^ - ■ ■ , vi) , with L the total number of single particle states, is a state of v unpaired bosons {v = v\) 
defined by 

K~ \v) = , n\ \v) = v\ \v) , 

v\ are referred as the Seniority quantum numbers. 

The total number of particles is N = 2M + v, with M the number of paired bosons. Each eigenstate l^) is completely 
defined by a set of M spectral parameters (pair energies) E m which are a particular solution of the Richardson's equations 

1=1 £(=im) = l 



and their eigenvalues are given by 



M 



1=1 



Hamiltonians constructed as general linear combinations of the integrals of motion H = Xa=i e i Riiv) have eigenvalues 
E = J2i=i £ i ( r i + ^1) where n is the eigenvalue of the integral of motion R\ 



B. Pairing Hamiltonians for bosons in confining traps 

As an application of the boson rational family, we will consider the problem of a boson system confined to a harmonic- 
oscillator trap and subject to boson pairing interactions l4lll . The pairing Hamiltonian with uniform couplings cannot describe 
the physics of a trapped boson system, for the following reason. Looking back at the commutators of the pair operators, Eq. 
d37i . we see that they are proportional to the degeneracy of the state 1 that appears inside the definition of the generator Kf 
(Eq. I8QI >. Thus, the matrix elements of the pairing Hamiltonian between states 1 and 1' will be proportional to y/rt\fl\i. In a 3D 
harmonic oscillator with 1 = (nlm), the shell degeneracy is fl\ ~ n 2 . On the other hand, the single-particle energies are £1 = n. 
Thus, the net effect would be the scattering of boson pairs to high-lying levels with greater probability than to low-lying levels, 
producing unphysical occupation numbers. This was precisely the behavior observed in a numerical solution of Richardson's 
equations J83l > for a system of 1000 bosons with an attractive pairing strength g l4lll . 

We can use the freedom we have in choosing the parameters 771 entering in the definition of the R\ operators to obtain 
a physically relevant exactly-solvable model. In order to cancel out the unphysical dependence of the pair-coupling matrix 
elements on the degeneracies, we choose the 771's so that 771 = (ei) 3 . The Hamiltonian, which is given by the linear combination 
of the new Ri's is 

H TB = 2 Y eiRi = C + E + E 5"' [ K \ K v - nxnv] , (85) 
1 1 1/1' 
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where 



l 1^1' 




1 



e\ + Ey + £\£y 



(86) 



The interaction in Eq. ( 186t has the nice feature that its two-body matrix elements decrease with the number of shells, as one 
would expect in general. It has the particular property that the interactions of the pair- and density-fluctuations are strictly the 
same but opposite in sign. Taking into account that e\ is proportional to n, the two-body matrix elements in Eq. ( I85l l cancel out 
the dependence on the degeneracies in the effective pair-coupling matrix elements. Thus, Hjb should be more appropriate than 
-Hbbcs when modelling a harmonically-confined boson system with a pairing-like interaction. 

The spectrum of Hjb can be obtained from the eigenvalues ?'i of the associated R\ operators as E = 2 J^i £ i r h with the end 
result being 



(Note that the first two terms of Eq. ( I87t exactly cancel the constant term C in Eq. j85ll). 

We solved Richardson's equations for Hjb for a system of M = 500 boson pairs and L = 50 oscillator shells. In this case, 
the occupation numbers display a reasonable physical pattern, with the occupancies decreasing monotonically with increasing 
single-boson energy 14111 . 

In the case of repulsive pairing a highly unexpected feature was found l4lll . For small values of g the system behaves as a 
normal BEC. At a critical value of the pairing strength g c a second-order quantum phase transition takes place. The new phase 
is characterized by a fragmentation of the condensate with the two lowest states macroscopically occupied, while the occupation 
of the other levels is negligible. 



The restriction of the bosonic RG models to two-levels comprises several well known quantum models. Among them we will 
discuss the Interacting Boson Model fl3ll . the LMG model |7] and the two Josephson-coupled BECs Hamiltonian fTill . Let us 
begin by defining the two integrals of motion from the most general XX Z RG models 



where the operators are defined in Eqs. d80t . From Eqs. (1881 we observe that the sum gives the total number of bosons which 
is a conserved quantity. We are then left with one independent quantum invariant, that we can take as the difference between the 
two integrals of motion to define the Hamiltonian 




(87) 



C. Exactly-solvable two-level boson Hamiltonians 



R a = Kl - X 12 [K+K- + K-K+] + 2Z 12 K z a K* 
R b = Kl + X 12 [K+K- + K-K+] - 2Z 12 K z a K z b 



(88) 



Hb2 =e(R b - R a ) = e (K* - K z a ) + 2e [X 12 (K+K- + K~K+) - 2Z 12 K* a K z b ] 



(89) 



Using Eq. d80t we rewrite Hb2 as 




where v = |X 12 and w = —sZ 12 are two arbitrary real numbers, and b^Ja^) creates a boson in level b(a) with an internal 
quantum number (3(a). As usual the bar in the internal labels means a time-reversed state, and £lb(a) is the degeneracy of the 
level. The constant term is C = | (fib — fl a — Zi 2 rtbCl a ). Using the parametrization of Eq. ( fT4l (with t\ = —r\ and t 2 = -q), 
4v 2 — w 2 = sg 2 e 2 , which can be positive, negative or zero depending upon the choice of parameters. 



The (unnormalized) eigenstates of Hb2 are given by 




(91) 
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with spectral parameters satisfying Bethe's equations 



ivd+Et + egd_(l + sEj) ^4 1 + s E t E n 

g~2r? + / I e - 1 + 2ge £ — '" =0, (92) 

where = d a ± d&, and 2d a ^ = v a {b) + ^a(b)/2- The corresponding eigenvalues can be constructed from the integrals of 
motion eigenvalues as 

w aaa ao ^ ^d-Ej + egd+{\ + sEj) 

Eb2 = id a db w - ed- - 2r\ ^ ^ . (93) 

1=1 1 

We will discuss next the application of the two-level RG bosonic models to three well-known quantum models. 



1. The Interacting Boson Model 

The Interacting Boson Model (IBM) has been a highly successful phenomenological model to describe the collective prop- 
erties of medium and heavy nuclei. The IBM captures the collective dynamics of nuclear systems by representing correlated 
pairs of nucleons with angular momentum L by ideal bosons with the same angular momentum. In its simplest version, known 
as IBM1, there is no distinction between protons and neutrons and only angular momentum L = 0(s) and L = 2(d) bosons 
are retained. The model has a U(6) group structure and three possible dynamical symmetry limits representing well defined 
nuclear phases: the U (5) symmetry for vibrational nuclei, the 0(6) symmetry for 7-unstable nuclei, and the SU (3) symmetry 
for axially deformed nuclei. In each of the three limits the Hamiltonian can be expressed in terms of the Casimir operators of 
the group decomposition chain. The three limits are then, exactly solvable with analytic expressions for the eigenstates. 

The transition from U (5) to 0(6) can be modelled by a boson paring Hamiltonian for the form 

1 — 2 

ffiBMi = x(n d - n s ) + N X ^ (d^d^ss + sVd^J , (94) 

fj,= -2 

where N is the total number of bosons and x is a parameter that interpolates between the linear Casimir operator of U(5), for 
x = 1, and the quadratic Casimir operator of 0(6), for x — 0. The Hamiltonian -H)bmi can be derived from Eq. d90i by making 
the following identifications: d = b, s = a, fid = 5> O s = 1, e = 2x, w = 0, and v = ■ 

The transition from the spherical vibrational phase (t/(5)) to the 7-unstable deformed phase was studied within the integrable 
model described by the Hamiltonian of Eq. (1941 14411 . It was found a second order quantum phase transition for a critical value 
of the control parameter x. In fact this is a unique point of second order phase transitions in the complete parameter space of the 
most general IBM Hamiltonian. The second-order character of the transition is related to quantum integrability. The ground state 
eigenvalue of iJiBMi is an analytic function of the control parameter x and, though level crossings are allowed due to quantum 
integrability, there are no level crossings in the low-energy spectrum. 



2. The Lipkin-Meshkov-Glick Model 

The LMG model has been extensively used for decades to simulate the phase transition from spherical to deformed shapes in 
finite nuclei. As introduced in section IVBl it is a schematic model describing the scattering of particle-hole pairs between two 
shells of different parity a. Though it was known for a long time that the model was quantum integrable, some analytic solutions 
were found only quite recently 114311 using the algebraic Bethe ansatz, after having mapped the model onto a Schwinger-boson 
representation. Here we will show that the LMG model is exactly solvable: After a Schwinger-boson representation of angular 
momentum operators, we will map the LMG model onto the two-boson integrable Hamiltonians of Eq. ( I90> . 

In the Schwinger mapping of the su(2) algebra the generators are expressed in terms of two bosons a and b as 

S+ = tfa=(S-y , S z = I (6+6 — o+a) = \{n b -n a ) , (95) 

with the constraint 

2S = tfb + a^a = n b + n a . (96) 
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Inserting Eq. d95l into Eq. ( 1741 we obtain a Schwinger-boson representation of the LMG Hamiltonian 

W e V + + + x W 

Hlmg = -ttS + -(n b - n a ) + — [b r b r aa + oJtfbb) + —n b n a . (97) 

We then recover the two-boson exactly solvable Hamiltonian of Eq. J90i with f2& = Sl a = 1, v = V/2N and w = W/N. In 
particular, the (unnormalized) eigenvectors and eigenvalues of -Hlmg are 

My xx &t& | 



= 1 



lx \E e + r) E e -r)J 



Elmg = E B2 -j, (99) 

where £32 is given in Eq. ( I93i . and the spectral parameters E( satisfy Eq. ( 192k Notice that, for each sector S, the number of 
bosons that enter in the expression for | v I')lmg is constrained to be 2S = n a + n b - 

We would like to emphasize that our expressions are compact forms valid for any arbitrary set of parameters in the Hamil- 
tonian. In particular, they are valid for the parameter range W 2 < V 2 , which correspond to the solutions not found in Ref. 
ll45ll . Since we have shown (using the weakly-interacting limit solutions) that this Bethe ansatz covers all possible eigenstates, it 
implies that the LMG is exactly solvable. 



3. Two coupled Bose-Einstein Condensates 



The Josephson effect, predicted more than forty years ago |43], describes pair tunneling between two superconductors through 
an insulating junction. An analogous effect can be realized with trapped ultracold bosonic gases in two different ways. In the 
first setup, two atomic condensates in the same atomic state are separated by a controllable potential barrier. In the second setup, 
atoms are condensed in two overlapping hyperfine states with an exchange mixing interaction. Both systems are described by 
the Hamiltonian 

H s = (cU+dtc) + (cWc + dWd) , (100) 

where c^and eft are left or right trap boson creation operators, or they create bosons in two different hyperfine states, depending 
upon the particular setup. Ej is the Josephson coupling exchanging bosons between the two states, and E c is the charging energy. 
This Hamiltonian has been recently exactly solved using the algebraic Bethe ansatz In fact, it is no more than another form 
of the two-level boson pairing Hamiltonian of Eq. d90t . It can be easily recast in the two-level form after performing the unitary 
canonical transformation 

c = (a - ib) , d = (a + ib) . (101) 
v 2 V 2 

Eliminating irrelevant constant terms, the Josephson Hamiltonian can be rewritten as 

Hi = §(n 6 - O - ^ (tftfaa + a)o)bb - 2n b n a ) , (102) 

which can be easily related to the LMG Hamiltonian, -Hlmg* by choosing e = 2Ej/N, V = —NE c /4: = —W. Therefore, the 
physics of the two coupled BECs is completely analogous to that of the LMG model. 



VII. MIXING REALIZATIONS AND REPRESENTATIONS OF THE GAUDIN FIELD OPERATORS 



We have already mentioned that our scheme for generating completely integrable models relies on finding different repre- 
sentations of the GGA. So far, we have simply concentrated on exactly solvable models where every single su(2) or su(l, 1) 
generator labeled by an index in the set T is equivalently represented. We still have the freedom to mix the representations of 
these generators and develop new exactly-solvable model Hamiltonians. The fact that the two su(2) realizations of Eqs. d76t and 
1771 are mutually commuting was also noted in 13^1 . They called them the spin and the charge realizations, respectively. This 
property implies that the elements of the spin and charge su(2) algebras act on orthogonal Hilbert spaces, allowing to define 
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an integrable RG model in each space separately. Thus, the following Hamiltonian 13^1 was proposed to study the interplay 
between pairing correlations and spin-exchange interactions 

ffch-s = $>7f + ^ £ (Xij(T.+T.- + t.-t+) + 2Z i3 r[T? + WijSi • S S ) , (103) 

where we can recognize the most general integrable pairing Hamiltonian, Eq. J45t . with an additional spin-exchange interaction, 
which is also integrable if the matrix W is derived from the rational Gaudin family as W\j = XL. The use of the rational 
family assures the conservation of the total spin quantum number as well as the third component of the total spin. The model 
Hamiltonian -f/ch-s can still accommodate a linear term in the spin variables representing a non-uniform magnetic field, or even 
more general XX Z Gaudin models can be implemented in the spin space at the cost of breaking the spin rotational symmetry. 

A numerical study of the interplay between pairing and exchange interactions in small metallic dots has been carried out 
in l50ll for systems with up to 30 levels. Even though the model used was fully integrable, the numerical results were mostly 
obtained by large scale diagonalization methods due to the complexity in solving the two sets of coupled Richardson's equations. 
The recently developed numerical techniques IBOII to solve efficiently these set of nonlinear equations may help to extend these 
studies to larger grains. 

The space orthogonality between the two su(2) fermionic realizations can also be exploited by defining different integrable 
models in the charge and spin sectors. For instance, it would be possible to mix a RG integrable Hamiltonian in the charge space 
with a Heisenberg or Haldane-Shastry model in the spin space. By mixing different realizations of the GGA, one can generate 
spin-fermion, spin-boson, or simply spin models with spins belonging to different irreducible representations. It turns out that 
this mixing-representations scheme might be useful to study decoherence and dynamic phenomena in open quantum systems 
where some degrees of freedom correspond to the system while the others (coupled in a particular way to the system) represent 
the thermal bath. In the following we illustrate these ideas starting with the generalized Dicke (GD) model Hamiltonian. 



A. Generalized Dicke models 



The GD models have been recently derived from the XX Z RG models by replacing one of the su(2) copies by a single boson 
satisfying the Heisenberg-Weyl algebra ll47ll . This procedure can be rigorously followed by expressing the su(2) generators in 
the Holstein-Primakoff representation as 



S+ = V2~S~o & f y 1 - ^ = (So") f . S o = -So + tfb 

where we have distinguished the particular copy of su(2) by the label 0, w (b) is the creation (annihilation) operator of a boson 
with \b, w] = 1, and So is the magnitude of the spin of this particular representation. Since the RG model is integrable for 
arbitrary spin values we can then analyze the limit So — ► 00, which implies the replacement 

5+ = v/2S^ bl Sq = v/2S7 b . 

We will skip here the derivation of the new class of integrable spin-boson models that can be found in j47ll and present the 
final form of the integrals of motion 

Ro = w6 t 6 + 2££ j S*/ + (tfSr + bS 



.1 

-2 



(104) 



Ri = wS? + \sZ — ~r [ s t s i + s i s t + 2 S?Sj] ~ V ( bf Sr + bS?) - 2e ; Sf . 

j(^i) £i 6j 

It can be readily verified that the set of operators of Eqs. J104l > are Hermitian, independent, and mutually commuting. 
Therefore, they constitute a new class of integrable spin-boson models. Though they have been derived from the trigonometric 
family of the RG models, the set of L operators Ri are identical to the rational family of RG models, except for the last two 
terms of su(2) which are essential for ensuring the commutation with the new bosonic integral of motion i? . 

Any function of these operators defines an integrable Hamiltonian. In particular, we can recognize Ro as a Dicke Hamiltonian 
describing the interactions of a multi-atom system with a single-mode radiation field. Moreover, a linear combination involving 
the whole set of integrals of motion, Eqs. ( 11041 . gives rise to more general integrable spin-boson Hamiltonians. 

Richardson's ansatz for the common eigenstates of the integrals of motion is 

o=l \ i=l ya 1 / 
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where L is the total number of su(2) spins and M the total spin and third component of the system. The M spectral parameters 
y a are particular solutions of the set of M of non-linear coupled Richardson's equations 

The corresponding expressions for the eigenvalues of Eqs. J104I are 

ro = E y <* ~ E S i £ ) > 



2 2 ^ £; - £j ^ V a - 2£; 

j(#i) J « y 



B. An exactly-solvable Kondo-like impurity model 

The effect of magnetic impurities in metals has been a subject of intense debate since the early 1930s. In 1964 J. Kondo El 
made significant progress by providing an explanation to the problem of the resistance minimum (as a function of temperature) 
in some metals such as Au. He recognized that the interaction of a single magnetic impurity with the conduction electrons is 
well represented by the s-d exchange Hamiltonian 

H 5d = J kk'( S '~ c kT c k'.l + S+c U c k'l + ^(Ck-r^'T - c kj. c k'|)) > ( 107 ) 

k,k' 

where S Z ,S ± are the spin operators representing the localized moment of magnitude S, while creates a conduction electron 
with momentum k and spin projection a. This type of interaction is characterized by terms in which the spin of the electron is 
flipped upon scattering with the impurity and are essential to understand the logarithmic contribution to the resistivity, and thus, 
its minimum. The simplest Hamiltonian representing the interaction of a localized moment with a band of itinerant electrons is 
the Kondo impurity model 

H K = J2^ n ^ + H sd ■ (108) 

k,(j 

It is important to emphasize that the case S = 1/2 , Jkk' = J/N, and linear (relativistic) dispersion has been exactly solved by 
Andrei and Weigmann using the Bethe ansatz 14911 . 

To find a new exactly-solvable single impurity Kondo-like model let us consider the RG constants of motion, where the 
localized spin of magnitude S is singled out 

R a = BS z + J2( x ai(S'S^ + S + Sr) + 2Z oi S z S^ , (109) 
j 

with electron spins given by 

S+ - c] TCji - (Sr)t , 5? = i(n JT - n n ) . (110) 
As mentioned above the commuting su{2) algebra 

rt = 4t c Il = > T ! = |(^JT + n U - 1) , (HI) 

is a gauge symmetry of Rq, thus one can write down the following exactly-solvable Gaudin model 

j 

= E e >JT + n n -1) + BS Z + J2 (ykk'(S~4 T c k'i + 5 + C [ iCk , T ) + 4 k ,S z (4 T c k' T - c ki c k'i)) , (112) 

j k,k' 
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where 

J fi B) = izZ^-^o^) , andc] CT = * £^4, (H3) 



J k 



To derive an exactly-solvable Kondo-like Hamiltonian one considers the particular case J^, = JXq^S^ and J kk , 



is the Fourier-transformed electron operator. 

To derive an exactly-solvable Kondo-like 
t/Zok<5kk' in i? 5 d, and adds the conduction band term 

H GK = ^e k n kff + BS Z + ,/]T ( x ok(S~cl ]Cu + S+cl^) + Z ok ^(c kT c kT - c ki c ki )) , (114) 

k,<T k 

since in this case r k commutes with R (Notice that the magnitude of the localized spin is not restricted to S = 1/2). Both, 
a minimum in the electrical resistivity, and the formation of a singlet resonance state characterize the Kondo physics. Clearly, 
since Hqk is translationally invariant the impurity contribution to the charge resistivity is zero. However, using this model one 
may address the fundamental issue of the formation of the singlet state, writing down a many -body state that captures the essence 
of the Kondo problem. 

The (unnormalized) iV-particles eigenstates of Hqk are given by 

i*) = n ( x tos + +j2 x ^ c h c *i ) i ps > ' ( n 5) 

1=1 \ k / 

where | FS) is the tensor product of the state | v) = | f kl ■ ■ ■ ^ kj • • • ) (of v = ^ k ^ k paired (i/ k = 2) fermions) with the remaining 
N — v fermions in a ferromagnetic state, and the lowest-weight spin state |0}g 

|FS> = H®n c kll >®l°>S: (116) 



N-u 

while E/s satisfy the Bethe equations 

B M 

— =d Z Qe + J2d k Z ]ie + J2 Zne > e=l,---,M, (117) 

k n (^)=l 

and the energy eigenvalues are given by (i\ = 0, 1, 2) 

E = Y. £ k^ + d ° [B + 2jJ2 z oe + 2jJ2 d * Z( >A ■ ( H 8) 

k \ t k / 

The case of zero magnetic field (B = 0) corresponds to the Gaudin magnet. 

Any set of parameters Xok, ^ok for which X^ k — = T (see Eq. dl2» leads to an integrable exactly-solvable model. 
Particularly, the case of a spin-isotropic exchange interaction, i.e. Xq^ = Zok = l/(?7o — ?7k), is of interest since the total spin 
is a good quantum number. It is easy to see that if the exchange coupling is antiferromagnetic (e.g., J > 0, r/o > ?7 k , and B 
smaller than a critical value) a singlet Kondo many-body state emerges in the problem. The emergence of this state is formally 
connected to BCS superconductivity, and the connection is established through Bethe's equations. In this way, for M > 1, we 
are rigorously connecting the Kondo resonance with the Cooper resonance problems. 



C. Exactly-solvable spin-boson models 

The study of spin-boson systems, i.e. a single spin of magnitude S linearly coupled to a thermal bath represented by a set of 
harmonic oscillators, is of particular interest in the theory of open quantum systems. These systems display important features 
of decoherence, that is the dynamical loss of quantum coherence because of the environment. In this section we will describe 
two interesting spin-boson models that are exactly-solvable. 

Let us start from the constant of motion Rq of Eq. dl09t and add to it the total magnetization symmetry 

H sbl = (B+u)S z + J2Xos(S-S^ + S + Sr) + 2j2z jS z S^+iuJ2s? . (119) 
j j j 
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Following a similar procedure to the one illustrated in section lVII Al we represent the su(2) spins Sj (in the limit Sb — > oo) as 



St = bt Sr = J2S^ b, , Sf = -S b + b]b 



j ' j 



j J 



and choose 



X *= 9 1 + it 



with the resulting Hamiltonian (up to irrelevant constants) 



(120) 



(121) 



where V = g^/2Sb, and B = B + u>, with B = B + V ^ £j . Notice that the localized spin S may have arbitrary magnitude 
S. This model is known to be trivially solvable. Indeed, it is a particular case of the Frohlich-like Hamiltonian describing a spin 
coupled to longitudinal optical phonons 



ffsbi = BS> + Y l ublb J +Y f <y J *S-bl + V i S + b i ), 



(122) 



which is diagonalized by first performing a unitary canonical mapping to new bosonic modes Oj 



where(A=^£jl^jl 2 ) 



yielding the Hamiltonian 
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H sbl = BS Z +J2 u}a l a i + A (S~ a \ +S+ ai ) 



(123) 



(124) 



(125) 



representing, in the new basis, optical phonons effectively interacting with a single spin through a single mode. 

A more elaborate spin-boson model can be realized using the other constants of motion. In this way the following model 
Hamiltonian results 



H sb2 = BSz + ^tiMbi-btbd + V^S-bl + S+bs) 



(126) 



m, with £j representing the parameters defining the model, Eq. ( I120K and x\ the coefficients in the linear 



where ty = 
combination of the integrals of motion. 
The eigenstates of H s t>2 are 



where the spectral parameters y a are the solutions of the Bethe equations 

1 + £ sv<* _ 2V ^ 



Do 



and Oo = 2S + 1. The corresponding eigenvalues are 

~2 



^sb2 = -7tV ^y a - 1J 2^ 



2 £ j 
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VIII. DIFFERENTIAL OPERATOR REALIZATIONS: SCHRODINGER-GAUDIN OPERATORS 

So far we concentrated on discrete representations which led to various exactly solvable lattice models. However, it is well- 
known that it is possible to realize representations of Lie algebras in the form of differential operators. As we will see in 
this section these representations will lead to models in the continuum. Basically we will generalize the work pioneered by 
Ushveridze and others 0|. 

Out of the many applications one can foresee we will concentrate on a single problem. The problem consists of mapping 
exactly-solvable lattice models to their equivalent in the continuum. Clearly, in general, the full spectrum of the lattice will be 
embedded in the spectrum of the continuum equivalent. For the sake of simplicity we will use the ® j su(2) representation. 

To illustrate the procedure we will use the BCS pairing lattice model of Eq. (I58> which, after realizing about its underlying 
algebraic structure, reads 

#bcs = 2 Y <**k + 9 E T k rv ■ (127) 

k kk' 

We wonder what the many-body problem in configuration space (a continuous manifold), to which it maps onto, is. To this end, 
we will consider the corresponding su(2) Lie algebra of differential operators in the tensor product representation 

T k = z k, 

r k = -z k d 2 Zk + 2S k 8 Zk , (128) 
^ = ~Sk + z k d Zk , 

(129) 

with Casimir operator = |(7k rj" + t^t£) + t^t^ = Sk(Sk + 1) and where it is assumed that the differential generators 
act on polynomials in the variables z^ of maximum order 2Sk- In other words, the order of the polynomials depends upon the 
value Sk of the spin irreducible representation. In this differential operator representation Eq. d!27t can be written 

#bcs = 2Eek(-Sk + Zk^ k )+^({zk})E(-2k5 2 k +2S k 9 Zk ), (130) 

k k 

ffecs = E + gi p({z k })J22H9 Xk +A] 2 + V({x k }), (131) 

k 



where E = -2 ^ k £ kSk, </?({^k}) = Lk z k, z k = x\, 

(Sk + |) , Ek^k 



A = 2i( v - 4/ + f- * , and 4V({x k }) = i(8 Xk A) - A 2 . (132) 
V £k gip{{z k \)J 

Notice that, in this language, -Hbcs represents a many-particle system in a gauge field subjected to a potential V. 
In the case Sk = 1/2, for all k's (no unpair single -particle states) 

|T>k-*k, U)k-1 (133) 

and the representation space includes polynomials of degree at most 1 in each variable. In the following we will simply con- 
centrate on this case which corresponds to degeneracy 17k = 2. It is straightforward to prove that the ansatz wave function 

* M ({«k})= Y! s k 1 s k2 ---Sk AJ , (134) 

ki,"' ,k M 

where the prime in the sum means that ki ^ k-2 7^ • • • 7^ kjvf , and 

s km = Zk "' , (135) 

is a solution of Eq. Jl 3 II . Interestingly, the function v I /M ({sk}) represents an elementary symmetric function of order M. It 
turns out that these functions are the equivalent of the Richardson's solutions in the continuum, i.e, 

M 

#* M ({s k }) = (E + Y E m )* M ({s k }) , (136) 

m— 1 

with the complex numbers E m satisfying Richardson's equations 

1 M 1 

^ 2e k - E m j— 1 E m - Ei 
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IX. CONCLUSIONS 



The benefits of having exact solutions to problems involving strongly interacting many-particle systems are difficult to over- 
state. New exactly(or quasi-exactly)-solvable models are always a unique tool to better understand physical phenomena char- 
acterized by non-linear and non-perturbative effects. Moreover, exactly-solvable models are excellent testing grounds for ap- 
proximations to the many-body problem. In the present work we have explored a generalized Gaudin algebra, whose invariants 
provide the generating function for integrable quantum Hamiltonians called XYZ Gaudin models. These quantum invariants 
can be simultaneously diagonalized using the Bethe ansatz. Different representations of the generators of the generalized Gaudin 
algebra realize many well-known physical Hamiltonians including the Bardeen-Cooper-Schrieffer, Suhl-Matthias-Walker, Inter- 
acting Boson Model of nuclei, Lipkin-Meshkov-Glick, several BEC models, generalized Dicke, spin-boson, a new Kondo-like, 
and many other models not yet exploited in the literature. In this way, we have identified the underlying algebraic structure, thus 
providing a unifying framework. An advantage of the Bethe ansatz for the Gaudin models is that the physical interpretation of 
their eigenfunctions is straightforward. The built-in correlation physics is so transparent that they could have well been chosen 
as exact variational states. 

An important question concerns the differences between mean-field approximations to the eigenfunctions of the Gaudin 
models and their exact solutions in the thermodynamic limit. Here, we briefly discussed issues related to the quantum critical 
behavior of these models. In particular, we analyzed the nature of the transition between the superconducting and Fermi liquid 
phases in the BCS model. We concluded that it is of the Kosterlitz-Thouless type independently of the space dimensionality of 
the lattice. 

A number of applications have been presented with the intention of illustrating the variety of physics problems described by 
microscopic Hamiltonians which belong to the class of XX Z Gaudin models. We have shown that the Lipkin-Meshkov-Glick 
model, widely used in nuclear physics and, more recently, in connection with Quantum Information Theory, is exactly-solvable. 
Our proof seems to complete the work initiated in Ref. l45ll to the whole parameter space (including the sector W 2 < V 2 ). 
The exact solvability of the two-level boson Hamiltonians given in lVI Cl comprises three important models, the LMG, the IBM 
(describing the transition from the i/(5) to the 0(6) dynamical symmetries), and the two Josephson-coupled BECs. All of them 
are, therefore, characterized by the same physics. 

We have also shown that the XY Z Gaudin equation, Eq. (J6J, results from the use of the Jacobi identities for the generators of 
the algebra, i.e. it is a property of the algebra. For the XX Z case, we have derived a family of antisymmetric solutions which 
includes the rational, trigonometric, hyperbolic, and Richadson's solution as special cases. Indeed, we have proved that the latter 
is a reparametrization of any of the other three. 

Solving efficiently the non-linear Bethe equations, which provide the necessary spectral parameters, is an important technical 
issue. In this regard, knowing the strong- and weak-coupling limits of those equations help to reduce the complexity of the task. 
In the strong-interacting limit this analysis was previously done in Ref. 15111 . In AppendixlAlwe have shown that the weakly- 
coupled limit solutions of the Bethe equations are given by the roots of Laguerre polynomials, by transforming those equations 
into a generalized Stieltjes equation. In this way, we proved that the Bethe ansatz encompasses all possible eigenstates, and does 
not provide spurious solutions for finite couplings. 



APPENDIX A: THE WEAKLY INTERACTING LIMIT 



The Bethe equations, Eq. J40l >. can be expressed in terms of the parametrization of Eq. dl4> as 



l-2gs(j: ieT dj+M-l 



A I 



In the limit g — > one recovers the non-interacting model, with the variables tf converging to the parameter values tj, depending 
on the corresponding distribution of the particles or spins. Therefore one can make the substitution 

h = ti l +gxt, (A2) 

where \t is the index of the parameter value to which tf converges. Then Eq. JA1I can be rewritten up to first order in g as 

2d;, v-^ 2 1 + 2qa\, 

— + E = - /. J » < A3 > 



with a$ e independent of the variables xf. 

M 



E d 3 1 Prr-+ E 1 -^^-s(d h +N ie -l), (A4) 
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and where N$ t is the number of variables t n that cluster around the parameter t- u . 

Eq. (IA3I can be transformed into a Generalized Stieltjes equation 15 211 . with the solutions given by: 

1 + st? 

Xi = T~T~n H, (A5) 

where the r; are the roots of the associated Laguerre polynomials L k N , with fc = 2d- }l — 1 and AT = N- u . Note that the resulting 
values for the variables are correct up to second order in g: 

h = t h +g .^ l) n. (A6) 
1 + 2ga k 

The variables r\ will be real for d- u > (with su (1, 1) realizations, typical for bosons), while for dj c < (with su(2) realizations, 
typical for fermions) the variables xe will come in complex conjugated pairs, except for one real value in case of an odd number 
of variables per cluster. 

Because a polynomial of order N has a unique set of TV roots, the weakly interacting limit establishes a one-to-one mapping 
between the non-interacting solutions (defined by the number of variables clustered around each parameter ij) and the solutions 
at finite values of g. 



ACKNOWLEDGMENTS 



This work was supported by the US DOE through Contract No. W-7405-ENG-36, by the Spanish DGI under grant BFM2003- 
053 16-C02-02 and by the Fund for Scientific Research - Flanders (Belgium). Fruitful discussions with C. D. Batista, V. G. Guer- 
guiev, I. Martin, S. Trugman and D. Van Neck are acknowledged. 



[1] T. Senthil, et al. Science 303, 1490 (2004). 

[2] M. Gaudin, La Fonction d'Onde de Bedie (Paris: Masson, 1983). 

[3] A. G. Ushveridze, Quasi-exactly solvable models in quantum mechanics (IOP, Bath, 1994). 

[4] CD. Batista and G. Ortiz, Phys. Rev. Lett. 85, 4755 (2000). 

[5] J. Bardeen, L. N. Cooper, and J. R. Schrieffer, Phys. Rev. 108, 1175 (1957). 

[6] H. Suhl, B. T. Matthias, and L. R. Walker, Phys. Rev. Lett. 3, 552 (1959). 

[7] H. Lipkin, N. Meshkov, and A. Glick, Nucl. Phys. 62, 188 (1965). 

[8] R. H. Dicke, Phys. Rev. 93, 99 (1954). 

[9] J. Dukelsky, C. Esebbag, and P. Schuck, Phys. Rev. Lett. 87, 066403 (2001). 
[10] L. Amico, A. Di L orenzo, A. Osterloh , Phys. Rev. Lett. 86, 5759 (2001). 

[11] R. W. Richardson, cond-mat/02035 1 2 

[12] C. D. Batista and G. Ortiz, Adv. in Phys. 53, 1 (2004); cond-mat/0207 lOfO 

[13] F. Iachello and A. Arima, The Interacting Boson Model (Cambridge University Press, Cambridge 1987). 

[14] S. Kohler, and F. Sols, Phys. Rev. Lett. 89, 60403 (2002) 

[15] R. Somma, G. Ortiz, H. Barnum, E. Knill, and L. Viola, quant-ph/0403035 

[16] R. Gilmore, Catastrophe Theory for Scientists and Engineers (Wiley-Interscience, New York, 1981). 
[17] M. Gaudin, J. Physique 37, 1087 (1976). 
[18] In the XYZ case they are given by 

cm _ y Q x QV — V QV c z — _ ST^ 7 ^ 2 
jer j6T jer 

[19] H. Kamerlingh Onnes, Comm. Leiden 119b, 120b, 122b, 124c (191 1); J. de Noble, Phys. Today 49 (9), 40 (1996). 
[20] See, for example, M. Tinkham, Introduction to Superconductivity (McGraw-Hill, New York, 1996). 
[21] A. Bohr, B. R. Mottelson, and D. Pines, Phys. Rev. 110., 936 (1958). 
[22] The XYZ Gaudin models will be analized somewhere else. 

[23] R. W. Richardson, Phys. Lett. 3, 277 (1963); Phys. Lett. 5, 82 (1963); Phys. Lett. 14, 325 (1965); J. Math. Phys. 6, 1034 (1965); Phys. 

Rev. 144, 874 (1966); Phys. Rev. 159, 792 (1967); R. W. Richardson and N. Sherman, Nucl. Phys. 52, 221 (1964). 
[24] M.C. Cambiaggio, A.M.F. Rivas and M. Saraceno, Nucl. Phys. A 424, 157 (1997). 
[25] J. von Delft, and D.C. Ralph, Phys. Rep. 345, 61 (2001). 
[26] J. M. Roman, G. Sierra, and J. Dukelsky, Phys. Rev. B 67, 064510 (2003). 
[27] D. Gobert, M. Schechter, U. Schollwock, J. von Delft, cond-mat/0404614 



24 



A THE WEAKLY INTERACTING LIMIT 



[28] R. W. Richardson, Phys. Rev. 141, 949 (1966) 

[29] G. Sierra, J. Dukelsky, G. G. Dussel, J. von Delft, and F. Braun, Phys. Rev. B 61, Rl 1890 (2000). 

[30] S. Rombouts, D. Van Neck, and J. Dukelsky, to be publish in Phys. Rev. C.lnucl-th/03 1 20701 

[31] R. W. Richardson, J. Math. Phys. 18, 1802 (1977). 

[32] J. M. Roman, G. Sierra, and J. Dukelsky, Nucl. Phys. B 634, 483 (2002). 

[33] H. R. Krishnamurthy and B. S. Shastry, Phys. Rev. Lett. 84, 4918 (2000). 

[34] J. Dukelsky, C. Esebbag, and S. Pittel, Phys. Rev. Lett. 88, 062501 (2002). 

[35] L. Amico, A. Di Lorenzo, and A. Osterloh, Nucl. Phys. B 614, 449 (2001). 

[36] F. Pan, J. P. Draayer and W. E. Ormand, Phys. Lett. B 422, 1 (1998). 

[37] S. Rombouts, D. Van Neck, K. Peirs & L. Pollet, Mod. Phys. Lett. A 17, 1899 (2002). 

[38] D. Van Neck, Y. Dewulf M. Waroquier, Phys. Rev. A 63, 062107 (2001). 

[39] C.K. Law, H. Pu, and N.P. Bigelow, Phys. Rev. Lett. 81, 5257(1998). 

[40] R.W. Richardson, J. Math. Phys. 9, 1327 (1968). 

[41] J. Dukelsky, and P. Schuck, Phys. Rev. lett. 86,4207 (2001). 

[42] J. Dukelsky, and S. Pittel, Phys. Rev. lett. 86,4791 (2001). 

[43] B. D. Josephson, Phys. Lett. 1, 251 (1962). 

[44] J.M. Arias, J. Dukelsky, and J.E. Garcfa-Ramos. Phys. Rev. Lett. 91, 162502 (2003). 
[45] F. Pan, and J.P. Draayer, Phys. Lett. B 451, 1 (1999). 

[46] H.-Q. Zhou, J. Links, R.H. McKenzie,and X.-W. Guan, J. Phys A 36, LI 13 (2003). 

[47] J. Dukelsky, G. G. Dussel, C. Esebbag, and S. Pittel, to be published by Phys. Rev. Lett. I cond-mat/0406190 1 
[48] J. Kondo, Prog. Theor. Phys.32, 37 (1964). 

[49] N. Andrei, Phys. Rev. Lett. 45, 379 (1980); P. B. Wiegmann, Sov. Phys. JETP Lett. 31, 392 (1980). 

[50] G. Falci, R. Fazio, and A. Mastellone, Phys. Rev. B 67, 132501 (2003). 

[51] E. A. Yuzbashyan, A. A. Baytin, B. L. Altshuler, Phys. Rev. B 68, 214509 (2003). 

[52] B. S. Shastry and A. Dhar, J. Phys. A: Math. Gen. 34, 6197 (2001). 



25 



